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ABSTRACT 

We present a method for reducing cosmological data to constraints on the amphtudes of modes of the 
dark energy density as a function of redshift. The modes are chosen so that 1) one of them has constant 
density and 2) the others are non-zero only if there is time-variation in the dark energy density and 3) 
the amplitude errors for the time- varying modes are uncorrelated with each other. We apply our method 
to various combinations of three- year WMAP data (Spergel et al., 2006), baryon acoustic oscillation data 
(Eisenstein et al., 2005), the Riess et al. (2004) 'Gold' supernova data set, and the Supernova Legacy 
Survey data set (Astier et al., 2005). We find no significant evidence for a time- varying dark energy 
density or for non-zero mean curvature. Although by some measure the limits on four of the time- 
varying mode amplitudes are quite tight, they are consistent with the expectation that the dark energy 
density does not vary on time scales shorter than a Hubble time. Since we do not expect detectable time 
variation in these modes, our results should be viewed as a systematic error test which the data have 
passed. We discuss a procedure to identify modes with maximal signal-to-noise ratio. 

Subject headings: cosmology: theory - cosmology: observation 



1. INTRODUCTION 

That the cosmological expansion rate is accelerating is 
well-established^. Two pressing questions, of deep impor- 
tance for fundamental physics, are 1) Is the acceleration 
due to corrections to the gravitational field equations or 
to some unknown matter component that we call dark 
energy? and 2) If dark energy, is it a cosmological con- 
stant or something with a time-varying density? Here we 
present a method for reducing cosmological measurements 
of distance as a function of redshift in a manner suited to 
answering this second question. 

Our method reduces cosmological data (anything that 
depends on the history of the dark energy density, Pxi^)) 
to constraints on a cosmological constant plus the am- 
plitudes of time- varying modes of Px{z). The modes are 
chosen to have uncorrelated amplitude errors and to be 
those that are best determined by the data. The two chief 
desirable properties of such a reduction are: 

1. The cosmological constraints can be expressed in 
a highly model-independent manner in terms of 
just a few numbers (plus associated functions of 
redshift). 

2. Consistency with a cosmological constant is 
straightforward to study visually from the graphed 
results. 

Many have previously considered diff'erent methods for 
parameterizing possible departures from a cosmological 
constant (Chevalher & Polarski, 2001; Weller & Albrecht, 
2002; Linder, 2003; Huterer & Starkman, 2003). Such pa- 
rameterizations, while not necessary for comparing the rel- 
ative merits of two dark energy models, facilitate interpre- 
tation of the data in a less model-dependent manner. Most 
of these are one or two-dimensional parameterizations of 

^See for example Shapiro & Turner (2005) and references therein. 



w(z). 

The most frequently used parameterizations are w(z) is 
constant and w{z) = wq + {1 — a) Wa where a = 1/(1 -I- z) 
is the scale factor (Chevallier & Polarski, 2001; Linder, 
2003). However, even allowing for non-zero Wa, these al- 
low only for departures from constant density in a highly 
restricted space among the space of all possible ways the 
density could vary. It is possible for experiments to be sen- 
sitive to time variation and yet have wq — —1 and Wa — 
to within the uncertainties. Additionally, current data do 
not constrain wq and Wa very well, and Simpson & Bridle 
(2006) find that the parameterization can lead to signif- 
icant biases in the inferred w{z) and other cosmological 
parameters. 

We have chosen to work with density because it has a 
more direct relation to the data. The data we consider 
depend on distance as a function of redshift. This func- 
tion can be calculated from Px{z) with one integral while 
calculating it from w{z) requires two integrals. 

Wang & Freese (2006), Daly & Djorgovski (2004) and 
Wang & Tegmark (2004) reconstruct Px{z) from distance 
and other cosmological data. Although these reconstruc- 
tions are useful for visual inspection, detailed interpreta- 
tion of results is obscured by the correlation of errors across 
redshift. To ameliorate this difficulty, Huterer & Cooray 
(2005) reduce data to weighted averages of Px{z) with un- 
correlated errors. This has advantages, but the drawback 
that to look for time-variation one must visually differen- 
tiate the data. 

Although we are working with Pxiz), we are not at- 
tempting a reconstruction of this function from the data. 
We identify the best-determined time-varying modes of 
Px{z) and then determine the probability distribution of 
their amplitudes, having marginalized over all the other 
parameters, including the amplitude of the constant mode. 
If any of these mode amplitudes is significantly non-zero, 
we have evidence for time- varying dark energy. By design, 
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the ampUtudcs of those errors are uncorrelated for ease of 
interpretation. 

Somewhat similar to our approach has been the study 
of eigenmodes of w{z) (Huterer & Starkman, 2003). The 
cigcnvahic spectra (Song & Knox, 2004; Knox ct al., 2005) 
suggest that w{z) may eventually need to be described 
with more than just two numbers, although for another 
view see Linder & Huterer (2005). 

Finally, we should mention that Wang & Tegmark (2005) 
reduce data to H(z) in rcdshift bins with the attractive 
property that the errors in H{z) are uncorrelated from 
bin to bin and only dependent on the supernovae in that 
bin. However, this reduction is not ideal for detecting a 
departure from a cosmological constant, especially given 
uncertainty in and 0^ which affect the conversion of 
H{z) to px{z). 

Our method is applicable to any measurements and any 
model space in which dark energy's sole influence comes 
through the history of its energy density, Px{z). We ap- 
ply it here to two different supernova data sets (the Gold 
(Riess et al., 2004) and SNLS (Astier et al., 2005)) with 
and without the baryon acoustic oscillation (BAO) dis- 
tance constraints from Eisenstein et al. (2005) and always 
with constraints from CMB data. The CMB data are im- 
portant for how it constrains the distance to last scattering 
and the matter density. 

Our method may be even more attractive as descriptions 
of the uncertainties in luminosity distances become more 
complex (due to the increased importance of systematic 
relative to statistical error). In this case the use of lumi- 
nosity distances as a final stage of data reduction becomes 
more cumbersome. A reduction to the constraints on a 
few mode amplitudes may then be significantly easier to 
use than the reduction to constraints on hundreds or thou- 
sands of luminosity distances. 

Identifying well-determined modes is useful not only for 
detection of time-variation, but also for ferreting out sys- 
tematic errors. With the best-determined modes, one can 
split the data up into subsets and check the corresponding 
mode amplitude estimates for consistency. Since the am- 
plitudes have small errors, they can survive the greatest 
amount of sub-sampling while remaining sufficiently well- 
measured that the consistency tests remain meaningful. 

We compare the data-determined limits on the mode 
amplitudes to theoretical upper limits derived from the 
criterion that the time-scale for dark energy density vari- 
ation is longer than a Hubble time. We find these Hubble 
time upper limits to be tighter than the data-determined 
upper limits for all modes. Thus, any detection of non- 
zero mode amplitude with current data would most likely 
be an indication of systematic error. Due to this theo- 
retical expectation of a null result we conclude that our 
method, as currently implemented on current data, is best 
though of as a powerful method for identifying otherwise 
undetected systematic errors. We discuss means of folding 
in theoretical expectations for the level of time variation so 
that we can identify modes that maximize signal-to-noise 
ratio, rather than just minimize noise. 

In our applications we make no assumption about the 
mean curvature. We allow it to vary, constrained only by 
the data, and plot the result together with 0,^ and the 
amplitude of the time- varying dark energy density modes. 



We do this to avoid what would be the highly unfortu- 
nate mistake of declaring detection of time variation when 
the data could just as well be explained by non-zero mean 
curvature. Linder (2005) has recently demonstrated the 
extent of the fix ^ wq — Wa degeneracy for future CMB 
-|- supernova data, and the degeneracy between and 
Wo for the SDSS BAO data is commented on in Eisenstein 
et al. (2005). We also see determination of the (possibly 
non-zero) mean curvature as another very interesting ap- 
plication of distance vs. rcdshift measurements (Bernstein, 
2005; Knox, 2006; Freivogel et al., 2005). 

In section II we describe the method in detail. In section 
HI we apply it to the SNLS + CMB data. Results from 
other combinations of supernova, CMB and BAO data are 
presented in Appendix B. In section IV we discuss our 
results and conclude. 

2. METHOD 

The goal of our analysis is to demonstrate a method for 
measuring non-constant dark energy with a combination of 
low-z distance measurements, CMB data, and BAO con- 
straints. In 2.1 we describe our parameterization of the 
cosmology. In 2.2 wc describe how the likelihood of these 
parameters is calculated given each of the data sets. In 
2.3 we describe our calculation of the eigenmodes. In 2.4 
we describe our use of the Monte Carlo Markov Chain 
method that we use to estimate the parameters and their 
uncertainties. 

2.1. Parametrization 

Our set of cosmological parameters is the matter density 
today, Urn in units of 1.88 x 10^^^ g/cm^, Qk = —kc?/Hl = 
1 — f2tot and the Ui that determine the history of the dark 
energy density Px{z). Given the basis functions, ei{z) (to 
be described in 2.3), the dimensionless parameters a, spec- 
ify Px{z) up to an overall constant which is the critical 
density today, pc- These cosmological parameters are sum- 
marized in Table 1. 

The supernova data sets also include some 'nuisance' pa- 
rameters required to model the data because of our inabil- 
ity to infer precisely the luminosity of each of the super- 
novae. These parameters arc described in the appropriate 
likelihood calculation subsections (2.2.3 and 2.2.4). 

2.2. Likelihood Calculation 

Our constraints come from combinations of WMAP3 
data (Spergel et al., 2006), a BAO constraint from Eisen- 
stein et al. (2005), and supernova data from Riess et al. 

Table 1 
Cosmological Parameters 



Matter density: lo^ = n^h^ h = lookm^ec/Mpc 
flk Curvature: fi^ = —k-^ = 1 — fltot 
Q.i Dark energy parameters: Px{z) = PcY^iCy.iei{z) 
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(2004) and Astier et al. (2005). Taking uniform priors on 
the parameters, the probability distributions of the pa- 
rameters given the data are proportional to the likelihood 
functions, L. Since we make no use of the probability 
amplitudes here, wc use likelihood and probability inter- 
changeably in what follows. 

2.2.1. CMB Likelihood 

CMB data are sensitive to a large number of cosmo- 
logical parameters, but only two of them are relevant for 

interpreting data that is sensitive to the expansion rate 
at lower rcdshifts, such as supernova data and BAO data. 
These are and the comoving angular-diameter distance 
to last scattering, = Duiz*) where 

^ , . 1 r rr f" dz' 

and z* is the redshift of the surface of last scattering, here 
taken as the peak of the visibility function. We therefore 
derive, from the CMB data, a likelihood distribution for 
the two-dimensional parameter space: lu„i,L)*- 

To calculate this likelihood function we use an MCMC 
chain calculated from the WMAP3 data set and available 
on the LAMBDA archive^. We use a chain which uses the 
simplest-case ACDM model with only WMAP data, which 
assumes zero curvature and a cosmological constant. How- 
ever, at fixed u)m and D* the CMB data are highly insen- 
sitive to departures from a cosmological constant or non- 
zero curvature. A chain that allowed for curvature and 
dark energy would not give a significantly different likeli- 
hood function for and as long as the dark energy 
and curvature remained sub-dominant at last scattering. 
The dominant effect of curvature or dark energy on the 
CMB is to change the projection of length scales on the 
last-scattering surface to observed angular scales because 
of the inffuence of k and H{z) on Dm{z*) in Eq. 1. The 
constraint from CMB data on Dm('Z*) thus captures the 
CMB information about f2fc and dark energy. 

Prior to the recent WMAP3 release we calculated the 
likelihood of ^Imh^ and from an MCMC chain calcu- 
lated in Chu et al. (2005) from WMAPl, Cosmic Back- 
ground Imager (CBI) and Arcmimite Cosmology Bolome- 
ter Array Receiver (ACBAR) data the 'WMAPext' data 
set used in Spergel et al. (2003). Our update to WMAP3 
does not have a significant impact on our dark energy re- 
sults. 

There is some weak information from the CMB about 
dark energy and curvature, beyond that in Um and D^, 
coming from the integrated Sachs Wolfe (ISW) effect (for a 
review of CMB physics see Hu & Dodelson, 2002). Ideally 
our CMB likelihood function would capture this informa- 
tion. For simplicity we ignore it, because the constraints 
are very weak and model-dependent. 

To calculate our CMB likelihood function, we begin by 
calculating the parameters D* and cOm at every chain step 
to create a new chain in these two parameters. By count- 
ing the number of chain elements in each bin of a regular 
grid in cOm and £>* we create a two-dimensional matrix 
that is an approximation to the probability distribution 
PcMsit^m, I)*) given by the WMAP3 data. Due to the 

■^Legacy Archive for Microwave Background Data Analysis: 
http: / /lambda.gsfc.nasa.gov 



properties of MCMC, the number of chain elements that 
fall within each bin will be proportional to the probability 
in that bin plus sample variance fluctuations. Since we 
have ignored the other chain parameters in this process, 
such as the optical depth to last scattering and the pri- 
mordial power spectrum spectral index, we have, in effect, 
marginalized over them. 

Before making use of this matrix, we attempt to re- 
duce the noise caused by the finite number of samples 
by implementing a low-pass filter. We do so by taking 
a Fast Fourier Transform of the two dimensional matrix 
above and eliminating the high frequency components. 
The inverse Fourier transform of this new matrix is a 
noise-reduced approximation to the probability distribu- 
tion. An appropriate frequency cut is chosen with two 
criteria: the difference between the original matrix and 
the noise-reduced matrix can be well-explained by Poisson 
noise, and the majority of the noise is eliminated (that is, 
neighboring bins do not vary dramatically). 

Due to the inaccuracy of both the original chain and 
this noise-reduced matrix for very low probabilities, we 
set to zero all bins which include fewer chain steps than 
some cutoff value, with this ciitoff value chosen so that 
fewer than 0.003 times the total number of chain elements 
are excluded when these bins are set to zero. If this zero- 
ing is not done, then the low probability regions of the 
noise-reduced matrix are dominated by ringing and in- 
clude negative values. The result of reducing the noise 
on an example probability distribution PcMBi'jJmiDt) is 
shown in Fig. 1. This noise-reducing algorithm allows us 
to have an approximation to the full chain that is closer to 
the chain output than a Gaussian approximation, has no 
obvious noise, and does not have bins that are so coarse 
that detail is lost. We evaluate Pcmb(wto, -D*) by bilinear 
interpolation over the smoothed grid. 

We have not included some other recent CMB results, 
most importantly the 2003 flight of Boomerang (Montroy 
et al., 2005; Piacentini et al., 2005; Jones et al., 2005) 
and the most recent results from CBI (Readhead et al., 
2004a, b). Including these would further tighten up the 
CMB constraints on LOm and Z?* (MacTavish et al., 2005; 
Sievers et al., 2005). Although, given the lack of change 
from WMAPext to WMAP3, we do not expect that these 
additional data would be significant for our application. 

2.2.2. BAO Likelihood 

The next step is to incorporate the BAO data from lumi- 
nous red galaxies in the SDSS survey. For this we use the 
parameter Dy given in Eisenstein et al. (2005) equation 

(2): 

Dy{z) ^ [Dl,{z)^Y'\ (2) 

Eisenstein et al. (2005) have compressed the full data 
set into a function of Dy in two separate ways. For this 
paper, we select the A parameterization, their equation 
(4), because it appears to give tighter constraints on the 

dark energy: 

A = Dv{Q:ih)^^^^^ = 0.469 ± 0.017(3.6%). (3) 

Though this BAO data is a very powerful constraint, it is 
unfortunate for our analysis that it was reduced to a con- 
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II = m — M = Slog 



10 



lOpc 



(7) 



For the Gold data set, the distance modulus n is esti- 
mated for each supernova. But since the absolute magni- 
tude M is unknown, we consider the /x for each supernova 
to be: 

IXiiSM) = f,f- SM, (8) 
where /if is given for each supernova in the Gold data set 
from Riess et al. (2004). The fif have been estimated by 
Riess et. al. from observations in order to correct the 
observed magnitude of each supernova for dust, as well 
as differences in magnitude between supernovae that are 
correlated with the shape of the luminosity versus time 
function. In doing this estimation to generate the Gold 
data set, a mean value of the absolute magnitude M has 
been assumed. The parameter 5M is the difference be- 
tween the true mean absolute magnitude of the supernovae 
and the estimated absolute magnitude. This parameter is 
marginalized over. To evaluate for a cosmology with a 
distance-redshift relation Dl{z) we calculate: 



X 



E 



(M.(<5M)-51ogio (^)y 



(9) 



Fig. 1. — Probability distribution of cOm and D» given CMB data 
before smoothing (top) and after smoothing (bottom). Note that 
the large-scale features are preserved, while the small-scale, sample- 
vaxiance-induced noise has been greatly reduced. 



straint at a single redshift. This reduction was intended to 
be valid for the case of a cosmological constant, and fairly 
robust against changes in constant w{z), but may intro- 
duce significant systematic errors for the cases we consider 
with much more freedom in the behavior of the dark en- 
ergy. It would be much more useful to have the data re- 
duced to distance constraints for multiple redshift bins, 
and we encourage such reductions from future analyses of 
BAO data. 

We write = — 21n£ for the BAO data set as: 
{A - 0.469)2 



Xbao — 



0.0172 



(4) 



2.2.3. Gold Likelihood 



The observed magnitude m of an object can be linked 
to cosmology by: 

'Dl{z) 



m = M + 51ogiQ 



lOpc 



(5) 



where the absohite magnitude M is defined as the hypo- 
thetical observed magnitude of an identical object at lOpc, 
and the luminosity distance Dl{z) is evaluated at the red- 
shift of the object in question. The luminosity distance is 
simply-related to the comoving angular-diameter distance 
defined in equation 1 by a factor of {1 + z): 



Dl{z) = {1 + z)Dm{z). 



(6) 



We can then define the distance modulus n as the por- 
tion of equation 5 that is dependent only upon the cos- 
mology as: 



This assumes Gaussian noise in the distance moduli 
quoted in the Gold data set, and the sum over i is over 
the various supernovae. 

2.2.4. SNLS Likelihood 

The SNLS data are reduced differently, providing us 
with constraints on the parameters used to calculate the 
effective apparent magnitude rather than just the effective 
apparent magnitude itself. For each supernova, the stretch 
factor s, color factor, c, and rest frame B band apparent 
magnitude are estimated from the light curves. We 
can then define a distance modulus /j. for each supernova 
which can be directly compared to the cosmology: 

lii{M,a,0)=m*B-M + a{si-l)-0Ci. (10) 
For our likelihood calctilation we take the best fit values 
of to|j, s, and c from Astier et al. (2005) for each super- 
nova, and then treat M, a, and f3 as global parameters 
free to vary just like the cosmological parameters. 

Since we do not have a covariance matrix for the param- 
eters m*g, s, and c, we simply fix these quantities to their 
best fit values for each supernova and use the uncertainty 
in /xb provided in Astier et al. (2005) to form 

^_^ (/x.(M,a,/3)-51ogio(^))' 

Here af^^ is the intrinsic dispersion of the absolute mag- 
nitudes. We take the value of (Tint = 0.15 for the nearby 
supernovae, and (Tint = 0.12 for the SNLS supernovae, the 
mean values reported in Astier et al. (2005). 

Residuals of the binned SNLS distance modulus data, 
after subtraction of the mean distance moduhis for the 
ACDM model, are shown in Fig. 3. The binning is purely 
for plotting purposes. For calculation of the likelihood 
function we use the unbinned data. 
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2.2.5. Final Likelihood 

We now have three independent data sets: CMB, BAO, 
and one of two supernova data sets. Since they are inde- 
pendent, the probabihties multiply: 



ln(£TOT) = ln(PcMB(Wm,-D*)) 



2 2 
XSN XbAO 



(12) 



2 2 

whore Xsw comes from cither the Gold or SNLS data set. 
Because the SNLS and Gold data sets make use of many of 
the same nearby supernovae, but use different algorithms 
for correcting the magnitudes, the data sets are correlated. 
Combining them appropriately would require more than a 
simple summing of the Fisher matrices and is an exercise 
we have not attempted. 

2.3. Eigenmode Calculation 

To select our modes we need the covariance matrix for 
the errors in the ai parameters for some initial basis. One 
could in principle calculate the covariance matrix from an 
MCMG chain, but since a chain run with the full set of ai 
parameters is highly non-Gaussian, the calculated covari- 
ance matrix is not a good approximation and does not give 
good eigenmodes. Instead we estimate it analytically from 
the Fisher matrix, and then use this estimate to define the 
modes. 

2.3.1. The Fisher Matrix 

First, we generate a Fisher matrix by expanding about 
a particular point. The point chosen is the mean point of a 
chain run with the aforementioned likelihood function, but 
using instead of ai. This Fisher matrix is calculated 
in two parts. For the supernova part we obtain: 



dxk 1 dxk 
dpi al dpj 



(13) 



where for SNLS pk = {u>rn,^k,M,a,/3,ao,ai,a2, ...,Q;ioo) 
and 

Xk = 51ogio (^^^) +M- a{sk - 1) + /3cfc (14) 

and for Gold pk = {c^m, ^k, ^M, uq, ai,a2, aioo) and 

'DL{zk) 



Xk = 51og 



10 



lOpc 



5M. 



(15) 



Since the BAO data arc written as a single constraint in 
A, which is a function of the chain parameters, calculating 
its Fisher matrix is even simpler: 



pBAO ^dAl_dA 
dpi (7^ dpj 



(16) 



Once this is calculated, we need to factor in the CMB 
data. Because performing numerical derivatives on our 
two-parameter probability matrix would be numerically 

unstable, we first calculate a covariance matrix from the 
chain calculated in section 2.2.1 which includes only iVm 
and Df,, then invert this covariance matrix to obtain a 
Fisher matrix, which we will call i^^^^. We convert this 
Fisher matrix to one in the parameters of interest as fol- 
lows: 



f^CMB 
■ ij 



F 



dp, 



(17) 



where v is the vector defined by w = [wmj^**]- Finally, 

F%^^ = F^^ + F^,^^ + F^,^^ . (18) 

This total Fisher matrix is then inverted to give us a co- 
variance matrix for the eigenmode calculation. The covari- 
ance matrix, once diagonalized, gives us our eigenmodes. 

2.3.2. Initial Basis and Diagonalization 

Our specific goals for the desired modes are that there 
be one and only one constant mode, and that the time- 
varying modes have amplitudes with uncorrelated errors. 
We can meet these two criteria by careful selection of the 
initial basis and diagonalization proccdiu'c. 

The initial basis must have the following properties: 

1. It must have one element which is a constant. 

2. The elements must be linearly independent. 

Note that the second property is satisified by any basis 
by definition. There are many possible bases that sat- 
isfy these criteria. We have chosen ours with the addi- 
tional (somewhat arbitrary) criteria that their dot prod- 
ucts, 'Y^a^ii'^ai^ji^a) — ^^iji where N is the dimension- 
ality of the space given by the discretization procedure. 
The normalization factor. N , keeps the shapes and ampli- 
tudes of the best-determined modes independent of N as 
the continuum limit is approached. For completeness, we 
describe the initial basis in Appendix A. 

Each basis element is defined based on its value at A'^ 
values of z between and 2. With the choice of uniform 
spacing in ln(0.01 -I- z) we find we are able to approach the 
continuum limit with a smaller value of N than with other 
possible choices, such as uniform spacing in z. The results 
in this paper were generated with = 50. Our basis is 
piecewise-linear in ln(0.01 -I- z). 

From the procedure in Appendix A we get N ai param- 
eters with which to describe px{z)- To obtain our eigen- 
modes we invert the Fisher matrix in the parameter space 
of these ai plus all our other parameters to obtain a co- 
variance matrix. Then we take the i > subspace of this 
covariance matrix. Diagonalizing this N — 1 hj N — 1 
covariance matrix gives us our non-constant eigenmodes, 
each of which has errors uncorrelated with the other eigen- 
modes. Including the constant mode completes the space. 
We call this new basis e|™(z), and its coefiicients a|™. 

We typically truncate this new space by only keeping 
the first few best-measured eigenmodes. This results in 
an MCMC chain which is both highly Gaussian and easy 
to interpret, since a significant departure from zero of any 
one of the first few af^ parameters is a signature of non- 
constant dark energy. In the following we omit the super- 
script "em" from the parameters to reduce notational 
clutter. 

Note that the error in the constant mode has not been 
decorrelated with the errors in the non-constant modes. 
We discuss this choice in Appendix A. 

2.4. Generating the Chain 

We use the Metropolis-Hastings algorithm to generate 
a Monte Carlo Markov Chain (Gamerman, 1997) as de- 
scribed in Christensen et al. (2001). Due to the low com- 
putational requirements we are able to run very long chains 
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with 6,010,000 elements. We then ignore the first 10,000 
steps and thin the chain by taking every 20*^* element, 
resulting in a chain of length 300,000. 

We can be sure that this thinned, 300,000-length chain 
has converged by comparing the parameter estimates be- 
tween different subsets of the chain. For example, if we 
take the subset of the first 5,000 elements of the thinned 
chain with 6 dark energy parameters for the SNLS 4- CMB 
data set combination, and compare those parameter esti- 
mates to those of the last 5,000 elements of the chain, those 
parameter estimates usually vary by less than cr/lO. De- 
pending on the selection of elements, occasionally one of 
the 10 chain parameters will vary by as much as cr/5. Since 
the chain has converged to an accuracy within about cr/10 
at the 90% confidence limit for each parameter after 5,000 
elements, and we run the chain to 300,000, convergence is 
not an issue. 

We use a generating function which is a multivariate 
normal distribution. Its covariance matrix is obtained 
from a pre-run. This pre-run is done in two iterations, 
each of length 100,000 with thinning by 20. The second 
iteration uses the covariance matrix from the first, with the 
covariance matrix from the second giving the covariance 
matrix used to run the final chain. This iterative process 
gives us a generating function that is a good Gaussian 
approximation to the full probability distribution. This 
match between generating function and posterior makes 
for an efficient exploration of the posterior. 

We set a (weak) prior constraints on ri,„ by bounding 
it to the interval < ftm < 1- 

3. RESULTS 

We start by examining one data set combination in de- 
tail: SNLS with CMB data. The first few eigenmodes 
are shown in figure 2. As is typical with such eigenmode 
decompositions, the frequency of oscillations tends to in- 
crease as the mode number increases. 

The interpretation of these modes is simplified by exam- 
ining them in the distance modulus space. The effect of 
varying a; on distance modulus is complicated by the fact 
that ai is correlated with LUm, ^k-, oiq, and the supernova 
parameter M. Thus in Fig. 3 we take, for example, the 
parameter ai in the center panel, and set it to its mean 
value plus cr(ai). We then examine the chain in a small 
region about this value to obtain the mean values for all 
of the other parameters in the chain and calculate 
for this point in the parameter space. Finally, we subtract 
the best-fit model for comparison. We also perform 
the same procedure at the mean value minus cr(ai) and 
then repeat for ao and 0^2 . 

Note that these modes have significant support at 2; < 
0.1 and even z < 0.02 and thus are affecting distances 
at redshifts as small as 0.02. This is due to the nearby 
sample that spans 0.015 < z < 0.125. The low-redshift 
support in these modes makes them potentially sensitive 
to large bulk flows (Shi, 1997; Shi & Turner, 1998; Zehavi 
et al., 1998) that systematically shift the redshifts of the 
observed supernovae from their Hubble flow values. Hui & 
Greene (2005) find surprising sensitivity of w determina- 
tions to peculiar velocity effects in forecasts for future ob- 
servations. Given the weight some of our best-determined 
modes place on low z data, peculiar velocities could poten- 



tially be contributing significant systematic errors to our 
results as well. We address this concern below. 

The whole spectrum of (square roots of) eigenvalues for 
these modes is shown in Fig. 4. The stars are the result 
of our Fisher matrix calculation and the triangles from 
MCMC. Including more than four modes in our MCMC 
calculation leads to gross degeneracies that do not occur 
in the Gaussian approximation; we have not been able to 
reliably calculate the spectrum beyond this first handful 
of modes. For these lowest modes we see the Gaussian 
approximation provides a description of the errors good 
to about 15%. For the Gold data (see Appendix B) the 
Gaussian approximation is even better. 

At first blush, the large number of low noise modes con- 
tradicts what we know about the eigenvalue spectrum of 
w(z). Knox et al. (2005) found that for a future space-based 
mission only a few w{z) modes could be reconstructed with 
errors smaller than 0.1. Here, for current data, we see 
eight time- varying px{z)/pc modes can be reconstructed 
with errors smaller than 0.1. Others have pointed out 
that Px{z) can be reconstructed with much smaller frac- 
tional errors than w{z) (Wang & Freese, 2006). It is well 
known that differentiating data makes them noisier and, 
roughly speaking, reconstructing w{z) requires one more 
differentiation of the data than a reconstruction of px{z). 

We only mean to point out here that there is this differ- 
ence, and that the reasons for it can be understood. We are 
not using this difference as an argument that one should 
consider only Px{z) and not w{z). Although w{z) has 
larger error bars than px{z)/pc, this is not because more 
information is lost in a reduction to w{z). There is a Px{z) 
for every w{z), obtainable by integration. What happens 
in integrating noisy modes olw{z) to get the corresponding 
Px{z) is that the rapid oscillations with larger departures 
from the fiducial value of -1 get averaged out, suppressing 
the rapid oscillations. A large amplitude (compared to 
fiducial value of -1) noise fluctuation on a w{z) mode will 
be a much smaller fluctuation on the corresponding Px{z) 
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Fig. 2. — The three best-measured eigenmodes for the SNLS -|- 
CMB data set. 
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Fig. 3. — The observable consequences of altering the mode coefficients. The points with error bars are the binned SNLS residuals after 
subtraction of the mean ACDM distance modulus. The curves are the residual mean distance moduli with different constraints on Oi with 
i = (left panel), i = 1 (center panel) and i = 2 (right panel). The constraints are that is held fixed at 1 cr above or 1 cr below its mean 
value (dashed curves). As the mode number is increased, the oscillation frequency increases. 



(compared to fiducial value of ^^pc). For discussion of the 
relative merits of parameterizing dark energy by Px{z) and 
w{z) see Linder (2004). 
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Fig. 4. — The one a errors for the first 40 of our 50 eigenmode 
amplitudes plotted for both the MCMC analysis and the Fisher 
matrix analysis for SNLS + CMB. The I — a errors from the Fisher 
matrix are the square roots of the eigenvalues of the inverse of the 
Fisher matrix associated with the respective modes. The modes 
are normalized to have length y/N where N is the number of basis 
elements, as described in 2.3.2. 



In figure 5 we plot our estimates of cto ~ 0-7, and 
some of the varying at parameters. We subtract 0.7 from 
ao for plotting convenience. Recall that f2A = Q^o when 
the time-varying modes have zero amplitude. The multi- 
ple data points for a given parameter are the results for 
different numbers of parameters held fixed. 

Note that the first three nonconstant mode amplitudes 
are consistent with zero. Like many others (Doran et al., 
2005; Daly & Djorgovski, 2005; Xia et al., 2005; Ichikawa & 
Takahashi, 2005; Sanchez et al., 2006; Nesseris & Perivolaropou- 
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Fig. 5. — Parameter estimates for the SNLS + CMB data set. For 
each parameter, we estimated the value multiple times, each time 
allowing a different number of dark energy parameters to vary. 



los, 2005) we are finding consistency of available supernova 
plus CMB data with a flat ACDM universe. 

Note also that the curvature constraints are robust to 
variation of these first few modes. They are sufhciently 
well-constrained that they do not lead to significant con- 
fusion with the curvature. 

We also see that the amplitudes of the nonconstant 
eigenmodes themselves do not change appreciably when we 
change how many mode amplitudes we vary. Within the 
Gaussian approximation this independence is expected, 
since the modes were chosen to have uncorrelated Fisher 
matrix amplitude errors. 

While we have shown that for the first few dark energy 
parameters our Gaussian approximation remains valid and 
the curvature constraint remains robust, we expect that if 
we add more and more dark energy parameters, eventually 
these facts will change. We expect this for two reasons. 
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First, the freedom to vary the more poorly constrainted 
parameters will allow us to move far enough away in pa- 
rameter space (from the point at which the Fisher matrix 
was evaluated) for the Gaussian approximation to break 
down. Second, even within the Gaussian approximation, 
the errors in the amplitudes of the dark energy modes are 
only constructed to be uncorrelated with each other; they 
are correlated with the errors in all of the other parame- 
ters, including the displayed flk and ao- Thus including 
more modes will weaken these constraints. In figure 6 we 
see the results of including up to mode number 6. 

SNLS + CMB 
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Fig. 6. — As we allow more freedom in px{z) by allowing more 
to vary, the errors in the inferred cosmological parameters grow. The 
leftmost estimate for each parameter comes from a chain allowing 
up to ag to vary. For ao — 0.7 this left-most estimate is off the chart. 

The breakdown of Gaussianity is very sudden. The mo- 
ment the chain includes the parameter ag, the errors on 
all other parameters increase dramatically. We believe the 
breakdown is driven by a correlation between the ag error 
and the ao error, which then allows for variations far from 
the fiducial model. With these large variations the mean 
values of our cosmological parameters end up describing 
a cosmology that is inconsistent with other experiments. 
For example, when we include the ae parameter, gets 
dragged up against the prior's upper bound of 1 so that 
0.89 + 0.10 - 0.27 at 95% confidence. 

Results from the other data set combinations are not 
dramatically different to those from SNLS plus CMB that 
we have considered here in detail. For completeness, we 
include them in Appendix B. 

We now turn to the concern that because our best- 
determined modes have strong support at very low red- 
shift they may be contaminated by (unmodeled) peculiar 
velocity- induced redshifts. We address it by constructing 
modes (as described in Appendix A) that only have vari- 
ation in the region 0.2 j z j 2.0. This region contains all 
the SNLS supernovae and none of the supernovae from the 
nearby sample. 

Our best-determined modes resulting from this proce- 
dure are shown in Fig. 7 and the corresponding parameter 



constraints are shown in Fig. 8. These mode amplitudes 
have larger error bars than the previous ones, but are safer 
from contamination by peculiar velocities. 

4. DISCUSSION 

Our modes have been defined based on which ones can 
be best determined by the data. We now turn to consid- 
eration of the expected signal contribution to these mode 
amplitudes. 

It is difficult to see how the dark energy density could 
vary on time scales faster than a Hubble time. Were the 
dark energy density to start diluting with the expansion 
as a" we would have dlnp/dt = nH; i.e. the time-scale 
for variation would be l/\nH\. Even were the dark en- 
ergy density to suddenly thermalize as massless radiation 
we would only have n — —4. Expressing n in terms of 
variation of the equation-of-state parameter away from its 
value for a cosmological constant {Sw = w— (—1)) we have 
n = SSw. Thus we roughly expect \dlnp/dt\ < H. 

Ignoring contributions to H from curvature and matter, 
we can write dhipx/dt/ H = dpx/ Px/da/a which can be 
conveniently calculated for our modes by finite difference. 
Note that the time-scale for variation is not a property 
of the mode itself, but a property of the mode and the 
amplitude of the mode. We evaluate it for each mode with 
an amplitude equal to its la uncertainty. We find they all 
have regions of redshift for which they exceed unity, as can 
be seen in the two panels of Fig. 9. 

Note that even though the far region modes in Fig. 7 
look smoother than those in Fig. 2, their amplitude con- 
straints are weaker with the result that the density vari- 
ation time scales, evaluated with an amplitude of Icr, are 
actually smaller than those for the full modes. 

We have also explored the idea of finding modes that 
maximize signal-to-noise ratio rather than the ones that 
merely minimize noise. We have attempted to do so using 
the technique of signal-to-noise eigenmodes Bond (1995); 
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Fig. 7. — The three best-measured eigenmodes for the SNLS 
-I- CMB data set with the restriction that the modes only vary at 
z > 0.2. 
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Fig. 8. — Parameter estimates for the SNLS + CMB data set 
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each parameter, we estimated the value multiple times, each time 
allowing a different number of dark energy parameters to vary. 



Bunn & Sugiyama (1995); Karhunen (1947), also known 
as the Karhunen-Loeve technique. 

The procedure in general is for data with both a noise 
covariance matrix and a signal covariance matrix; i.e., data 
that is modeled as 

A, = Si+ n, (19) 
where the signal contribution has covariance matrix 



((.,-(.,))(.,■-(.,-))) 



(20) 



and similarly for the noise, Nij = {riinj). 

The procedure starts with a whitening of the data so 
that all modes have errors that are uncorrelated with unit 
variance. The signal matrix in this new space is N'^/^SN'^^"^. 
Eigenmodes of this 'signal-to-noise matrix' have eigenval- 
ues that are equal to the signal-to-noise ratio which is ex- 
actly what we want. 

For our application we view the density in redshift bins 
as the 'data'. The trick is identifying a useful form of the 
signal matrix. To date we have not found a useful form. 
For the signal matrices which we have tried, the result has 
been that though the modes themselves become smoother, 
the error estimates increase, such that the density varia- 
tion time scale, evaluated with an amplitude of la, does 
not decrease substantially. 

5. CONCLUSIONS 

The single most important question about the dark en- 
ergy is whether it is a cosmological constant. We have 
presented a method to test this hypothesis by looking for 
the deviations from constant density that are easiest to 
detect. The best-determined time- varying modes of Pxiz) 
are found by use of a Gaussian approximation to the pa- 
rameter likelihood, but the estimates of the amplitudes of 
these modes do not use the Gaussian approximation. We 
have applied the method to supernova, CMB and BAO 
data. 
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Fig. 9. — (Apx / Px) / {^ci/o,) for our modes with amplitudes at 
their 68% confidence upper limits. This is roughly equal to the ratio 
of the natural time-scale for variation (l/H) to the mode variation 
time-scale (din p/dt)~^ . We see that the mode variation time-scales 
are longer than the natural time scale. Top panel: modes varying 
over all z. Bottom panel: modes varying only in the SNLS region. 



We found that SNLS -I- CMB data are capable of con- 
straining the amplitudes of three time-varying modes of 
Px{z) to better than 5%. None of them are significantly 
different from zero. We found that if only these well- 
constrained modes are allowed to vary then the curvature 
is also not significantly different from zero. We found that 
the best-determined mode amplitudes have highly inde- 
pendent errors, which within the Gaussian approximation 
occurs by design. Wc found that they are sensitive to den- 
sity variations at very low (z ^ 0.02) redshifts. We noted 
that these redshifts are sufhciently small to raise concerns 
about the effects of bulk flows. We pointed out that the 
best-determined modes are useful not only for detection of 
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time variation, but also for ferreting out systematic errors 
by studying how the ampHtude estimates fluctuate with 
varying data subsample. 

We pointed out that, by quite general considerations of 
the time-scales expected for dark energy density variation, 
at the achieved sensitivity levels we do not expect a de- 
tectable signal. Because we strongly expect a null result, 
our application can be viewed as a consistency test which 
has the potential for exposing systematic errors. This sit- 
uation is analogous to that with the (thcorctially unex- 
pected) B-modes which are a useful diagnostic in cosmic 
shear observations (e.g. Crittenden et al., 2002). As the 
statistical constraining power of the data improves we will 
begin to probe the regime where signals can be expected. 
We also suggested that with the adoption of an appro- 
priate signal covariance matrix we could use the signal- 
to-noise eigenmode technique to maximize signal-to-noise 
ratio rather than simply minimizing noise. 

In Appendix B our results show that the BAO data add 
significant additional constraining power to our cosmolog- 
ical parameter estimates. Unfortunately, the convenient 
reduction to a single number, though valuable for other 
applications where the dark energy is assumed to have a 
constant equation-of-state parameter, is not clearly appli- 
cable to our problem. It would be useful to have a reduc- 
tion of the BAO data to constraints on the distances to 
several redshifts, or as a single constraint on some linear 
combination of redshifts. 

We thank D. Huterer for a useful conversation that stim- 
ulated this work as well as A. Albrccht and M. Hudson for 
useful conversations. This work was supported in part by 
NASA grant NAG5-11098 and NSF grant 0307961. 
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Here we describe our construction of the initial basis 
that we use as the first step in creating the time-varying 
modes with uncorrelated amphtude errors. 

The first basis vector is a constant, which before nor- 
malization looks like (1, 1, 1, ... , 1). The next ^ basis 
vectors take the form (0, 0, 1, -1, 0, 0). As long as 
none of the non-zero values in this second set of basis vec- 
tors overlap with one another, they will be orthogonal to 
one another, and each is clearly orthogonal to the constant 
vector. 

To fill out the space further with a third set of basis vec- 
tors we notice that each of the second set of basis vectors is 
orthogonal to the constant vector, but only in steps of two. 
A new basis vector that maintains the same value in steps 
of two will be orthogonal to each previous basis vector: (1, 
1, -1, -1, 0, 0). To find the next one we consider that 
this basis vector is orthogonal to (1, 1, 1, 1, 0, 0), so 
with the previous argument all the basis vectors made so 
far will clearly be orthogonal to (1, 1, 1, 1, -2, -2, 0, 0). 
The next basis vector will then include 6 I's and two -3's, 
and so on. Careful counting reveals that there are — 1 
of this third set of basis vectors, which will complete our 
basis if we have an even number of elements. If we have 
an odd number of elements we can add a final basis vec- 
tor that takes the form (1, 1, 1, 1, -(n-1)), since every 
element of the second set of basis vectors contains a zero 
in the last position. 

Below is an example 9x9 basis (with vectors as rows), 
before normalization. 



1 




1 




-3 
1 



1 




-1 




-3 
1 



(Al) 



Each of these basis vectors defines an e°{z). To do this, 
we connect each element in the vector with a particular 
position in z. The value of e^(z) between these positions 
is given by linear interpolation. These basis vectors e^(z) 
are used to define Px{z) as in equation (A2): 



N 



(A2) 



i=0 



It is worth noting that while we demanded that the non- 
constant eigenmodcs were also zero-mean, there is an ar- 
bitrariness to the meaning of zero mean. It depends on 
the definition of the dot product, which is, in turn, basis- 
dependent. For example, the average value of a mode sam- 
pled in ln(z) is different from the average value of the same 
mode sampled in z. 

Given this arbitrariness, it is worth considering relaxing 
the zero-mean condition. With that condition abandoned, 
it becomes possible to use the new freedom to decorrelate 
the time- varying mode amplitude errors with the constant 
mode amplitude error. In principle one could do so in 
two steps: first add the right amount of constant to each 



of our time-varying modes to decorrelate it with the con- 
stant. Then re-diagonalize in the time-varying subspace. 
In practice we have been unable to circumvent numerical 
instabilities that occur when we attempt this procedure. 

In order to suppress variation of the near region, we 
chose to modify our pre-diagonalization basis. We select 
a basis where the first mode is a constant. The second is 
constant over the near region and the far region, with the 
two differring in sign, with the values of the near and far 
regions chosen to enforce orthogonality with the constant 
mode. The next group of modes varies only in the near 
region, with the far region set to zero. Its structure is 
the same as that described earlier for the varying basis 
elements. The final group of modes again has the same 
structure, but varies only in the far region. An example 
9x9 basis is shown below (before normalization). 
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The location of the split is chosen to coincide with the 

desired location before which no variation will occur. Once 
this basis is chosen, we carry out the diagonalization pro- 
cedure as described in 2.3.2, with the exception that we 
choose the subset of the covariance matrix that corre- 
sponds to those basis elements which only vary in the far 
region. 

APPENDIX B 

OTHER DATA COMBINATIONS 

Here we show results from application to other data 
combinations and comment on the difficulty of including 
the; BAO data. The results from the various other data 
set combinations we calculated are shown in figures 10, 11 
and 12. 

Adding in the BAO data, we find that we gain a fair 
amount of extra constraining power for both curvature and 
our dark energy parameters. For the SNLS data adding in 
the BAO data gives us the ability to measure one more 
dark energy parameter. Interpretation of these results 
though is severely hampered by the manner in which we 
incorporated the BAO data. As described in the text we 
utilized the Eisenstein et al. (2005) reduction of their data 
to a constraint on a combination of Dm{z) and H{z) at 
z = 0.35 and the matter density (their A parameter). Al- 
though they tested the validity of this data compression for 
constant w models, we are not sure how valid it is for the 
more general dark-energy model space that we consider. 
We include the BAO data only to show the potential sta- 
tistical power of that data for an analysis such as ours. 

Artifacts of the compression to A could be seen in the 
modes had we not employed the following trick: when cal- 
culating the Fisher matrix in order to get the modes, we 
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Reduction of Cosmological Data for the Detection of Time-varying Dark Energy Density 
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Fig. 10.- 
CMB data. 



Eigenmodes (left panel), Eigenvalue spectrum (center panel) and parameter estimates (right panel) for BAO plus SNLS plus 



held the H{z) in A fixed. If we were to include dH{z) / dui 
as a contribution to dA/dai then there would be a large 
spike in the modes at z = 0.35. This change to the eigen- 
modes makes no measurable difference to the parameter 
estimates. 
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Fig. 11. — As above, but for the Gold supernova data set combined with the CMB. 



_ ttg allowed to vary 
_ Oq, C(j allowed to vary 
_ a„, a,, Oj allowed to vary 
Up to Dig allowed to vary. 



1 to vary 
to vary. 



a„-0.7 a, 




Gold + BAG + CMB 



o 
u 

b 

I 



^ Gaussian approximation 
A Calculated from MCMC 



10 

mode number 



_ oLq allowed to vary 
_ Oq, Oj allowed to vary 
_ tto' ^i' allowed to vary 
_ Up to Og allowed to vary- 



-0.7 



Fig. 12. — As above, but for the Gold supernova data set combined vt^ith BAO plus CMB. 



